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I. INTRODUCTION 



Statistics is everywhere in Cosmology, today more than ever: cosmological data sets are 
getting ever larger, data from different experiments can be compared and combined; as the 
statistical error bars shrink, the effect of systematics need to be described, quantified and 
accounted for. As the data sets improve the parameter space we want to explore also grows. 
In the last 5 years there were more than 370 papers with "statistic-" in the title! 
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There are many excellent books on statistics: however when I started using statistics in 
cosmology I could not find all the information I needed in the same place. So here I have 
tried to put together a "starter kit" of statistical tools useful for cosmology. This will not 
be a rigorous introduction with theorems, proofs etc. The goal of these lectures is to a) he a. 
practical manual: to give you enough knowledge to be able to understand cosmological data 
analysis and/or to find out more by yourself and b) give you a "bag of tricks" (hopefully) 
useful for your future work. 

Many useful applications are left as exercises (often with hints) and are thus proposed in 
the main text. 

I will start by introducing probability and statistics from a Cosmologist point of view. 
Then I will continue with the description of random fields (ubiquitous in Cosmology), fol- 
lowed by an introduction to Monte Carlo methods including Monte-Carlo error estimates 
and Monte Carlo Markov Chains. I will conclude with the Fisher matrix technique, useful 
for quickly forecasting the performance of future experiments. 



11. PROBABILITIES 



A. What's probability: Bayesian vs Prequentist 

Probability can be interpreted as a frequency 

r-- (1) 

where n stands for the successes and N for the total number of trials. 

Or it can be interpreted as a lack of information: if I knew everything, I know that an 
event is surely going to happen, then = 1, if I know it is not going to happen then V — 
but in other cases I can use my judgment and or information from frequencies to estimate V. 
The world is divided in Prequentists and Bayesians. In general, Cosmologists are Bayesians 
and High Energy Physicists arc Prequentists. 

For Prequentists events are just frequencies of occurrence: probabilities arc only defined 
as the quantities obtained in the limit when the number of independent trials tends to 
infinity. 

Bayesians interpret probabilities as the degree of belief in a hypothesis: they use judg- 
ment, prior information, probability theory etc... 
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As we do cosmology we will be Bayesian. 
B. Dealing with probabilities 

In probability theory, probability distributions are fundamental concepts. They are used 
to calculate confidence intervals, for modeling purposes etc. We first need to introduce the 
concept of random variable in statistics (and in Cosmology). Depending on the problem at 
hand, the random variable may be the face of a dice, the number of galaxies in a volume 5V 
of the Universe, the CMB temperature in a given pixel of a CMB map, the measured value 
of the power spectrum P{k) etc. The probability that x (your random variable) can take a 
specific value is V{x) where V denotes the probability distribution. 

The properties of V are: 

1. V{x) is a non negative, real number for all real values of x. 

2. V{x) is normalized so that ^ / dxV{x) — 1 

3. For mutually exclusive events xi and X2, V{xi + X2) = V{xi) + V{x2) the probability 
of Xi or X2 to happen is the sum of the individual probabilities. V{xi + X2) is also 
written as V{xiUx2) or V{xi.OR.X2). 

4. In general: 

V{a,b)^V{a)V{b\a) ; V{b,a) ^V{b)V{a\b) (2) 

The probability of a and b to happen is the probability of a times the conditional prob- 
ability of b given a. Here we can also make the (apparently tautological) identification 
V{a, b) = V{b, a). For independent events then V{a, b) = V{a)V{b). 



Exercise: "Will it be sunny tomorrow?" answer in the frequentist way and in the 
Bayesian way^ . 

Exercise: Produce some examples for rule (iv) above. 

^ for discrete distribution / — s- ^ 

^ These lectures were given in the Canary Islands, in other locations answer may differ... 
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While Frequentists only consider distributions of events, Bayesians consider hj^otheses 
as "events", giving us Bayes theorem: 

where H stands for hypothesis (generally the set of parameters specifying your model, al- 
though many cosmologists now also consider model themselves) and D stands for data. 
V{H\D) is called the posterior distribution. V{H) is called the prior and V{D\H) is 
called likelihood. 

Note that this is nothing but equation [2] with the apparently tautological identity 
V{a, b) = V{b, a) and with substitutions: b — > H and a — > D. 
Despite its simplicity Eq. [3] is a really important equation!!! 

The usual points of heated discussion follow: "How do you chose V{H)7" , "Does the 
choice affects your final results?" (yes, in general it will). "Isn't this then a bit subjective?" 



Exercise: Consider a positive definite quantity (like for example the tensor to scalar 
ratio r or the optical depth to the last scattering surface r). What prior should one use? a 
flat prior in the variable? or a logaritmic prior (i.e. fiat prior in the log of the quantity)? for 
example CMB analysis may use a fiat prior in Inr, and in Z = exp(— 2r). How is this related 
to using a flat pror in r or in r? It will be useful to consider the following: effectively we 
are comparing V{x) with V{f{x)), where / denotes a function of x. For example x is r and 
fix) is exp(-2r). Recall that: V{f) = V{x{f)) " 
appears here to conserve probabilities. 



The Jacobian of the transformation 



Exercise: Compare Fig 21 of Spergel et al (2007) (l26l) with figure 13 of Spergel et al 
2003 (1251 ) . Consider the WMAP only contours. Clearly the 2007 paper uses more data than 
the 2003 paper, so why is that the constraints look worst? If you suspect the prior you are 
correct! Find out which prior has changed, and why it makes such a difference. 

Exercise: Under which conditions the choice of prior does not matter? (hint: compare 
the WMAP papers of 2003 and 2007 for the fiat LCDM case). 
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C. Moments and cumulants 



Moments and cumulants are used to characterize the probabihty distribution. In the 
language of probability distribution averages are defined as follows: 



These can then be related to "expectation values" (see later). For now let us just introduce 
the moments: /i^ = {x'^) and, of special interest, the central moments: /i^ — {{^ ~ (^))"^)- 

Exercise: show that fio — 1 and that the average {x) — fii. Also show that 112 — 



Here, H2 is the variance, /is is called the skewness, /i4 is related to the kurtosis. If you 
deal with the statistical nature of initial conditions (i.e. primordial non-Gaussianity) or 
non-linear evolution of Gaussian initial conditions, you will encounter these quantities again 
(and again..). 

Up to the skewness, central moments and cumulants coincide. For higher-order terms 
things become more complicated. To keep things a simple as possible let's just consider 
the Gaussian distribution (sec below) as reference. While moments of order higher than 3 
are non-zero for both Gaussian and non-Gaussian distribution, the cumulants of higher 
orders are zero for a Gaussian distribution. In fact, for a Gaussian distribution all moments 
of order higher than 2 are specified by /ii and 112- Or, in other words, the mean and the 
variance completely specify a Gaussian distribution. This is not the case for a non-Gaussan 
distribution. For non-Gaussian distribution, the relation between central moments and 
cumulants k for the first 6 orders is reported below. 




(4) 




At5 = «5 + IOK3K2 



A*2 = K2 



1^2, = «3 








(5) 
(6) 
(7) 
(8) 
(9) 
(10) 
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D. Useful trick: the generating function 



The generating function allows one, among other things, to compute quickly moments 
and cumulants of a distribution. Define the generating function as 

Z{k) — {exp{ikx)) — j dxexp{ikx)V{x) (11) 

Which may sound familiar as it is a sort of Fourier transform... Note that this can be written 
as an infinite series (by expanding the exponential) giving (exercise) 



m = E (12) 



n=0 ^■ 



So far nothing special, but now the neat trick is that the moments are obtained as: 

fln^i-n^Zik)\k=o (13) 

and the cumulants arc obtained by doing the same operation on InZ. While this seems 
just a neat trick for now it will be very useful shortly. 



E. Two useful distributions 

Two distributions are widely used in Cosmology: these are the Poisson distribution and 
the Gaussian distribution. 



1. The Poisson distribution 

The Poisson distribution describes an independent point process: photon noise, radioac- 
tive decay, galaxy distribution for very few galaxies, point sources .... It is an example of 
a discrete probability distribution. For cosmological applications it is useful to think of a 
Poisson process as follows. Consider a random process (for example a random distribution of 
galaxies in space) of average density p. Divide the space in infinitesimal cells, of volume 6V 
so small that their occupation can only be or 1 and the probability of having more than one 
object per cell is 0. Then the probability of having one object in a given cell is Vi — pSV and 
the probabihty of getting no object in the cell is therefore Vo — i — pSV. Thus for one cell 
the generating function is Z{k) — 7^^ exp(iA;n) = 1 + p5V{exp{ik) — 1) and for a volume 
V with V/SV cells, we have Z{k) = (1 + pSV{exp{ik) - 1))^/^^ ~ exp[pl/(exp(iA;) - 1)]. 
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With the substitution pV — > X we obtain 
Z{k) = exp[A(exp(i/c) — 1)] = X]^o -^"/'"■' ^^pI^-^) 6^P(^^^)- Thus the Poission probability 
distribution we recover is: 

P„ = — exp[-A] (14) 
nl 



Exercise: show that for the Poisson distribution (n) = A and that cr^ = A. 



2. The Gaussian distribution 

The Gaussian distribution is extremely useful because of the "Central Limit theorem". 
The Central Limit theorem states that the sum of many independent and identically dis- 
tributed random variables will be approximately Gaussianly distributed. The conditions for 
this to happen are quite mild: the variance of the distribution one starts off with has to be 
finite. The proof is remarkably simple. Let's take n events with probability distributions 
V{xi) and < Xj >= for simplicity, and let Y be their sum. What is V{Y)1 The generating 
function for Y is the product of the generating functions for the Xi. 



zy{k) = y: 



m=0 



ml 



-/i 



(15) 



for n — > oo then Zy(A;) — > exp[— < >]. By recalling the definition of generating 
function (eq. [TT]we can see that the probability distribution which generated this Z is 



V{Y) 



1 



exp 



1 



2 < x^ > 



(16) 



that is a Gaussian! 



Exercise: Verify that higher order cumulants are zero for the Gaussian distribution. 



Exercise: Show that the Central limit theorem holds for the Poisson distribution. 



Beyond the Central Limit theorem, the Gaussian distribution is very important in cos- 
mology as we believe that the initial conditions, the primordial perturbations generated from 
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inflation, had a distribution very very close to Gaussian. (Although it is crucial to test this 
experiment ally. ) 

We should also remember that thanks to the Central Limit theorem, when we estimate 
parameters in cosmology in many cases we approximate our data as having a Gaussian dis- 
tribution, even if we know that each data point is NOT drawn from a Gaussian distribution. 
The Central Limit theorem simplifies our lives every day... 

There are exceptions though. Let us for example consider N independent data points 
drawn from a Cauchy distribution: V{x) = [iia^l + [{x — x)/a]'^]^^. This is a proper 
probability distribution as it integrates to unity, but moments diverge. One can show that 
the numerical mean of a finite number of observations is finite but the " population mean" 
(the one defined through the integral of equation (jl]) with f{x) = x) is not. Note also that 
the scatter in the average of data points drawn from this distribution is the same as the 
scatter in 1 point: the scatter never diminishes regardless of the sample size.... 

III. MODELING OF DATA AND STATISTICAL INFERENCE 

To illustrate this let us follow the example from (28). If you have an urn with A^ red 
balls and M blue balls and you draw from the urn, probability theory can tell you what 
the chances are of you to pick a red ball given that you has so far drawn m blue and n red 
ones... However in practice what you want to do is to use probability to tell you what is the 
distribution of the balls in the urn having made a few drawn from it! 

In other words, if you knew everything about the Universe, probability theory could tell 
you what the probabilities are to get a given outcome for an observation. However, especially 
in cosmology, you want to make few observations and draw conclusions about the Universe! 
With the added complication that experiments in Cosmology are not quite like experiments 
in the lab: you can't poke the Universe and see how it reacts, and in many cases you can't 
repeat the observation, and you can only see a small part of the Universe! Keeping this 
caveat in mind let's push ahead. 

Given a set of observations often you want to fit a model to the data, where the model 
is described by a set of parameters a. Sometimes the model is physically motivated (say 
CMB angular power spectra etc.) or a convenient function (e.g. initial studies of large scale 
structure were fitting galaxies correlation functions with power laws). Then you want to 
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define a merit function, tliat measures tlie agreement between tlie data and tlie model: by 
adjusting tlie parameters to maximize the agreement one obtains the best fit parameters. 
Of course, because of measurement errors, there will be errors associated to the parameter 
determination. To be useful a fitting procedure should provide a) best fit parameters b) 
error estimates on the parameters c) possibly a statistical measure of the goodness of fit. 
When a) suggests that the model is a bad description of the data, then a) and b) make no 
sense. 

Remember at this point Bayes theorem: while you may want to ask: "What is the 
probability that a particular set of parameters is correct?", what you can ask to a "figure 
of merit" is "Given a set of parameters, what is the probability that that this data set 
could have occurred?". This is the likelihood. You may want to estimate parameters by 
maximizing the likelihood and somehow identify the likelihood (probability of the data given 
the parameters) with the likelihood of the model parameters. 

A. Chisquare, goodness of fit and confidence regions 



Following Numerical recipes ((|24|). Chapter 15) it is easier to introduce model fitting and 
parameter estimation using the least-squares example. Let's say that Di are our data points 
and y{xi\d) a model with parameters a . For example if the model is a straight line then a 
denotes the slope and intercept of the line. 

The least squares is given by: 

x' = Y.<D,-y{x,\a)]^ (17) 

i 

and you can show that the minimum variance weights are Wi = l/crf- 



Exercise: if the points are correlated how does this equation change? 

Best fit value parameters are the parameters that minimize the x^- Note that a numerical 
exploration can be avoided: by solving dx^/dai = you can find the best fit parameters. 
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1. Goodness of fit 



In particular, if the measurement errors are Gaussianly distributed, and (as in this ex- 
ample) the model is a linear function of the parameters, then the probability distribution 
of for different values of at the minimum is the distribution for u = n — m degrees 
of freedom (where m is the number of parameters and n is the number of data points. 
The probability that the observed even for a correct model is less than a value is 
P(x^ < x^jiy) = V{v/2,x^/2) = r(z//2,x^/2) where F stands for the incomplete Gamma 
function. Its complement, Q = 1 — P(z//2, x^/2) is the probability that the observed 
exceed by chance even for a correct model. See numerical recipes ^4) chapters 6.3 and 
15.2 for more details. It is common that the chi-square distribution holds even for models 
that are non linear in the parameters and even in more general cases (see an example later). 

The computed probability Q gives a quantitative measure of the goodness of fit when 
evaluated at the best fit parameters (i.e. at Xmin)- If <5 is a very small probability then 

a) the model is wrong and can be rejected 

b ) the errors are really larger than stated or 

c) the measurement errors were not Gaussianly distributed. 

If you know the actual error distribution you may want to Monte Carlo simulate synthetic 
data stets, subject them to your actual fitting procedure, and determine both the probability 
distribution of your statistic and the accuracy with which model parameters are recovered 
by the fit (see section on Monte Carlo Methods). 

On the other hand Q may be too large, if it is too near 1 then also something's up: 

a) errors may have been overestimated 

b) the data are correlated and correlations were ignored in the fit. 

c) In principle it may be that the distribution you are dealing with is more compact 
than a Gaussian distribution, but this is almost never the case. So make sure you 
exclude cases a) and b) before you invest a lot of time in exploring option c). 

Postscript: the "Chi-by eye" rule is that the minimum x^ should be roughly equal to 
the number of data-number of parameters (giving rise to the widespread use of the so-called 
reduced chisquare). Can you -possibly rigorously- justify this statement? 
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2. Confidence region 

Rather than presenting the full probability distribution of errors it is useful to present 
confidence limits or confidence regions: a region in the m-dimensional parameter space (m 
being the number of parameters), that contain a certain percentage of the total probability 
distribution. Obviously you want a suitably compact region around the best fit value. It 
is customary to choose 68.3%, 95.4%, 99.7%... Ellipsoidal regions have connections with 
the normal (Gaussian) distribution but in general things may be very different... A natural 
choice for the shape of confidence intervals is given by constant boundaries. For the 
observed data set the value of parameters do minimize the x^; denoted by Xmin- 
perturb a away from Sq the increase. From the properties of the distribution it 

is possible to show that there is a well defined relation between confidence intervals, formal 
standard errors, and Ax^- We In tabl^we report the Ax^ for the conventionals 1,2, and 
3 — 0" as a function of the number of parameters for the joint confidence levels. 

In general, let's spell out the following prescription. If fi is the number of fitted parameters 
for which you want to plot the join confidence region and p is the confidence limit desired, 
find the Ax^ such that the probability of a chi-square variable with fi degrees of freedom 
being less than Ax^ is p. For general values of p this is given by Q described above (for the 
standard 1,2,3— o" see table above). 

P.S. Frequentists use ^ lot. 

B. Likelihoods 

One can be more sophysticated than x^, if V{D) (D is data) is known. Remember 
from the Bayes theorem (eq|3]) the probability of the data given the model (Hypothesis) is 
the likelihood. If we set 'P(-D) = 1 (after all, you got the data) and ignore the prior, by 
maximizing the likelihood we find the most likely Hypothesis, or , often, the most likely 
parameters of a given model. 

Note that we have ignored V{D) and the prior so in general this technique does not give 
you a goodness of fit and not an absolute probability of the model, only relative probabilities. 
Frequentists rely on analyses where a goodness of fit can be established. 

In many cases (thanks to the central limit theorem) the likelihood can be well approxi- 
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mated by a multi-variate Gaussian: 



1 



cxp --Y.^D-y),C-\D-y), 



(18) 



(27r)"/2|(ieiC|V2 



where = {iPi — yi){Dj — i/j)) is the covariance matrix. 
Exercise: when are hkehhood analyses and analyses the same? 



1. Confidence levels for likelihood 

For Bayesian statistics, confidence regions are found as regions R in model space such 
that Jj^V{d\D)da is, say, 0.68 for 68% confidence level and 0.95 for 95% confidence. Note 
that this encloses the prior information. To report results independently of the prior the 
likelihood ratio is used. In this case compare the likelihood at a particular point in model 
space jC{a) with the value of the maximum likelihood Cmax- Then a model is said acceptable 



Then the threshold should be cahbrated by calculating the distribution of the hkelihood 

ratio in the case where a particular model is the true model. There are some cases however 
when the value of the threshold is the corresponding confidence limit for a with m degrees 
of freedom, for m the number of parameters. 

Exercise: in what cases?^ 



C. Marginalization, combining different experiments 

Of all the model parameters some of them may be uninteresting. Typical examples 
of nuisance parameters are calibration factors, galaxy bias parameter etc, but also it may 

^ Solution: The data must have Gaussian errors, the model must depend linearly on the parameters, the 
gradients of the model with respect to the parameters are not degenerate and the parameters do not affect 
the covariance. 



if 




(19) 
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be that we are interested on constraints on only one cosmological parameter at the time 
rather than on the joint constraints on 2 or more parameters simultaneously. One then 
marginalizes over the uninteresting parameters by integrating the posterior distribution: 

P{ai..aj\D)^ J daj+i,...damP{a\D) (20) 

if there arc in total m parameters and we are interested in j of them (j < m). Note that if 
you have two independent experiments, the combined likelihood of the two experiments is 
just the product of the two likelihoods, (of course if the two experiments are non independent 
then one would have to include their covariance) . In many cases one of the two experiments 
can be used as a prior. A word of caution is on order here. We can always combine 
independent experiments by multiplying their likelihoods, and if the experiments are good 
and sound and the model used is a good and complete description of the data all is well. 
However it is always important to: a) think about the priors one is using and to quantify 
their effects, b) make sure that results from independent experiments are consistent: by 
multiplying likelihood from inconsistent experiments you can always get some sort of results 
but it does not mean that the result actually makes sense.... 

Sometimes you may be interested in placing a prior on the uninteresting parameters 
before marginalization. The prior may come from a previous measurement or from your 
"belief. 

Typical examples of this are: marginalization over calibration uncertainty, over point 
sources amplitude or over beam errors for CMB studies. For example for marginalization 
over, say, point source amplitude, it is useful to know of the following trick for Gaussian 
likelihoods: 

P{a,..a^_,\DU -^-^ei-'^^^^^^^^^^ (21) 
•' {2tt) 2 ||6||2 



1 

^ ^ cxp 



v27r cr^ 



HA-A) 



2 

repeated indices are summed over and ||C|| denotes the determinant. Here, A is the am- 
plitude of, say, a point source contribution P to the Ce angular power spectrum, A is the 
m — th parameter which we want to marginalize over with a Gaussian prior with variance 
(T^ around A. The trick is to recognize that this integral can be written as: 

1 



P{ai..am-i\D) = Coexp 



-^Ci-2C2A + CsA^ 



dA (22) 
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(where Co.. .3 denote constants) and that this kind of integral is evaluated by using the 
substitution A — > A — C2/C3 giving something oc exp[— l/2(Ci — CI/C3)]. 
It is left as an exercise to write the constants explicitly. 

D. An example 

Let's say you wa„t to constrain cosmology by studymg clusters uumber couuts as a 

function of redshift. Here we follow the paper of Cash 1979 (3). The observation of a 
discrete number of clusters is a Poisson process, the probability of which is given by the 
product 

r = Ul,[erexp{-e.)/n,\] (23) 

where is the number of clusters observed in the i — th experimental bin and Cj is the 
expected number in that bin in a given model: Cj = I{x)6xi with i being the proportional 
to the probability distribution. Here 6xi can represent an interval in clusters mass and/or 
redshift. Note: this is a product of Poisson distributions, thus one is assuming that these 
are independent processes. Clusters may be clustered, so when can this be used? 

For unbinned data (or for small bins so that bins have only and 1 counts) we define the 
quantity: 

N 

C = -21nP = 2(E-^lnJi) (24) 

i=l 

where E is the total expected number of clusters in a given model. The quantity AC between 
two models with different parameters has a distribution! (so all that was said in the 
section applies, even though we started from a highly non-Gaussian distribution.) 

IV. DESCRIPTION OF RANDOM FIELDS 

Let's take a break from probabilities and consider a slightly different issue. In comparing 
the results of theoretical calculations with the observed Universe, it would be meaningless 
to hope to be able to describe with the theory the properties of a particular patch, i.e. to 
predict the density contrast of the matter 6{x) = 6p{x)/p at any specific point x. Instead, 
it is possible to predict the average statistical properties of the mass distribution ^. In 

^ A very similar approach is taken in statistical mechanics. 
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addition, we consider that the Universe we hve in is a random reahzation of all the possible 
Universes that could have been a realization of the true underlying model (which is known 
only to Mother Nature). All the possible realizations of this true underlying Universe make 
up the ensamble. In statistical inference one may sometime want to try to estimate how 
different our particular realization of the Universe could be from the true underlying one. 
Thinking back at the example of the urn with colored balls, it would be like considering 
that the particular urn from which we are drawing the balls is only one possible realization 
of the true underlying distribution of urns. For example, say that the true distribution has 
a 50-50 split in red and blue balls but that the urn can have only an odd number of balls. 
Clearly the exact 50-50 split cannot be realized in one particular urn but it can be realized 
in the ensamble... 



Following the cosmological principle (e.g. Peebles 1980 (23y), models of the Universe have 
to be homogeneous on the average, therefore, in widely separated regions of the Universe 
(i.e. independent), the density field must have the same statistical properties. 

A crucial assumption of standard cosmology is that the part of the Universe that we can 
observe is a fair sample of the whole. This is closely related to the cosmological principle since 
it implies that the statistics like the correlation functions have to be considered as averages 
over the ensemble. But the peculiarity in cosmology is that we have just one Universe, 
which is just one realization from the ensemble (quite fictitious one: it is the ensemble of all 
possible Universes). The fair sample hypothesis states that samples from well separated part 
of the Universe are independent realizations of the same physical process, and that, in the 
observable part of the Universe, there are enough independent samples to be representative 
of the statistical ensemble. The hypothesis of ergodicity follows: averaging over many 
realizations is equivalent to averaging over a large (enough) volume. The cosmological field 
we are interested in, in a given volume, is taken as a realization of the statistical process and, 
for the hypothesis of ergodicity, averaging over many realizations is equivalent to averaging 
over a large volume. 

Theories can just predict the statistical properties of 6{x) which, for the cosmological 
principle, must be a homogeneous and isotropic random field, and our observable Universe 
is a random realization from the ensemble. 

In cosmology the scalar field 6{x) is enough to specify the initial fluctuations field, and 
-we ultimately hope- also the present day distribution of galaxies and matter. Here lies one 
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of the big challenges of modern cosmology, (see e.g. ??) 

A fundamental problem in the analysis of the cosmic structures, is to find the appropriate 
tools to provide information on the distribution of the density fluctuations, on their initial 

conditions and subsequent evolution. Here we concentrate on power spectra and correlation 
functions. 

A. Gaussian random fields 

Gaussian random fields are crucially important in cosmology, for different reasons: first of 
all it is possible to describe their statistical properties analytically, but also there are strong 
theoretical motivations, namely inflation, to assume that the primordial fluctuations that 
gave rise to the present-day cosmological structures, follow a Gaussian distribution. Without 
resorting to inflation, for the central limit theorem, Gaussianity results from a superposition 
of a large number of random processes. 

The distribution of density fluctuations 6 deflned as^ S — Sp/p cannot be exactly Gaussian 
because the held has to satisfy the constraint 5 > — 1, however if the amplitude of the 
fluctuations is small enough, this can be a good approximation. This seems indeed to 
be the case: by looking at the CMB anisotropies we can probe fluctuations when their 
statistical distribution should have been close to its primordial one; possible deviations from 
Gaussianity of the primordial density fleld are small. 

If 5 is a Gaussian random fleld with average 0, its probability distribution is given by: 



where 5 is a vector made by the Si, denotes the inverse of the correlation matrix which 
elements are = (SiSj). 

An important property of Gaussian random flelds is that the Fourier transform of a 
Gaussian fleld is still Gaussian. The phases of the Fourier modes are random and the 
real and imaginary part of the coefficients have Gaussian distribution and are mutually 
independent. 

Let us denote the real and imaginary part of 5k by ReSk and ImS^ respectively. Their 

5 Note that (6) = 




VDetC"^ 




(25) 




exp 
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joint probability distribution is the bivariate Gaussian: 

P{Re6k, Im6k)dRe6kdIm6k = ^-exp 

where al is the variance in Re6]^ and ImSk and for isotropy it depends only on the magnitude 
of k. Equation fl2Bl) can be re- written in terms of the amplitude |5k| and the phase (p^. 

Pi\Sk\,(j)k)d\Sk\d(j)k = -^exp 

that is \Sk\ follows a Rayleigh distribution. 

From this follows that the probability that the amplitude is above a certain threshold X 

is: 

/oo 
— exp 
./X (J 

Which is an exponential distribution. 

The fact that the phases of a Gaussian field are random, implies that the two point 
correlation function (or the power spectrum) completely specifies the field. 



ReSl + Im6l 
2^1 



dRe6]idIm6]i 



(26) 



\Sk\d\Sk\d(t)k 



(27) 



2ai 



|(5k|(i|(5i( 



exp 



X 



(28) 



P.S. If your advisor now asks you to generate a Gaussian random field you know how to 
do it. (It you are not familiar with Fourier transforms see next section) 



The observed fluctuation field however is not Gaussian. The observed galaxy distribution 
is highly non-Gaussian principally due to gravitational instability. To completely specify a 
non-Gaussian distribution higher order correlation functions are needed^; conversely devi- 
ations from Gaussian behavior can be characterized by the higher-order statistics of the 
distribution. 

B. Basic tools 

The Fourier transform of the (fractional) overdensity field 6 is defined as: 

6i; = aJ d^r5{f) exp[-i^ ■ r\ (29) 



For "non pathological" distributions. For a discussion see e.g. (jlTI). 
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with inverse 

5(r)^B J d^k5j;ex.p[ik ■ r\ (30) 
and the Dirac delta is then given by 

6^{k) = BaJ d^r exp[±ik ■ r\ (31) 

Here I chose the convention A — 1, B — l/(27r)^, but always beware of the FT conventions. 
The two point correlation function (or correlation function) is defined as: 

C{x) = {5{r)5{f+ f )) = y < > exp[i^ • f] exp[i^ • (f + x)]d^kd^k' (32) 

because of isotropy ^(|a;|) (only a function of the distance not orientation). Note that in 
some cases when isotropy is broken one may want to keep the orientation information (see 
e.g. redshift space distortions, which affect clustering only along the line-of sight ). 
The definition of the power spectrum P{k) follows : 

< 5^5p >= {2nfP{k)S''{k + k') (33) 

again for isotropy P{k) depends only on the modulus of the k- vector, although in special 
cases where isotropy is broken one may want to keep the direction information. 
Since 5{r) is real, we have that 5t — S_^, so 

< 5^5*g >= (27r)3 J d^xCix) exp[-i^ • f]5'^(^ - k') (34) 



The power spectrum and the correlation function are Fourier transform pairs: 

e(^) = ^ / P{k)eM^k-f]d'k (35) 
P{k) = J ^{x) ex.p[-ik ■ x]d^x (36) 

At this stage the same amount of information is enclosed in P{k) as in ^{x). 
From here the variance is 

^2 =< S^x) >= e(0) = TTT^ / P{k)d^k (37) 



(27r)3 . 
or better 

(7^ = A\k)dlnk where A\k) = J^^^Pi^) (38) 
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and the quantity A^(/c) is independent form the FT convention used. 
Now the question is: on what scale is this variance defined? 

Answer: in practice one needs to use filters: the density field is convolved with a filter 
(smoothing) function. There are two typical choices: 

/ - ^2Tif/m% ^^P[- V2^V-Rg] Gaussian ^ a = exp[-A;2i?y2] (39) 

/ = (4;^®(^/^r) TopHat ^ /fe = ^^A_[sin(A;i?T) - kRTCOs{kRT)] (40) 
roue; hly Rt ~ \/^Rg- 

Remember: Convolution in real space is a multiplication in Fourier space; Mul- 
tiplication in real space is a convolution in Fourier space. 



Exercise: consider a multi-variate Gaussian distribution: 

where Cij —< 5i5j > is the covariance. Show that is 5i are Fourier modes then Cij is 
diagonal. This is an ideal case, of course but this is teUing us that for Gaussian fields the 
different k modes are independent! which is always a nice feature. 

Another question for you: if you start off with a Gaussian distribution (say from Inflation) 
and then leave this Gaussian field S to evolve under gravity, will it remain Gaussian forever? 
Hint: think about present-time Universe, and think about the dark matter density at, say, 
the center of a big galaxy and in a large void. 



1. The importance of the Power spectrum 

The structure of the Universe on large scales is largely dominated by the force of gravity 
(which we think we know well) and no too much by complex mechanisms (baryonic physics, 
galaxy formation etc.)- or at least that's the hope... Theory (see lectures on infiation) give 
us a prediction for the primordial power spectrum: 

PW^a(A)" (42) 
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n - the spectral index is often taken to be a constant and the power spectrum is a power 
law power spectrum. However there are theoretical motivations to generalize this to 



as a sort of Taylor expansion of n{k) around the pivot point /cq. dn/d\nk is called the 
running of the spectral index. 

Note that different authors often use different choices of feg (sometimes the same author 
in the same paper uses different choices...) so things may get confused.... so let's report 
explicitly the conversions: 



Exercise: Prove the equations above. 

Exercise:Show that given the above definition of the running of the spectral index, 
n{k) — n{ko) + dn/ din k\n{k / ko) . 

It can be shown that as long as linear theory applies -and only gravity is at play- S « 1, 
different Fourier modes evolve independently and the Gaussian field remains Gaussian. In 
addition, P{k) changes only in amplitude and not in shape except in the radiation to matter 
dominated era and when there are baryon-photon interactions and baryons-dark matter 
interactions (see Lectures on CMB) . In detail, this is described by hnear perturbation growth 
and by the "transfer function". 

C. Examples of real world issues 

Say that now you go and try to measure a P{k) from a realistic galaxy catalog. What are 
the real world effects you may find? We have mentioned before redshift space distortions. 
Here we concentrate on other effects that are more general (and not so specific to large-scale 
structure analysis). 




(43) 




(44) 



22 



1. Discrete Fourier transform 



In the real world when you go and take the FT of your survey or even of your simulation 
box you will be using something like a fast Fourier transform code (FFT) which is a discrete 
Fourier transform. 

If your box has side of size L, even if S{r) in the box is continuous, Sk will be discrete. 
The k-modes sampled will be given by 

k = (^^) {i,j,k) where = ^ (45) 

The discrete Fourier transform is obtained by placing the S{x) on a lattice of grid 
points with spacing L/N. Then: 

= ^^E^M-ik-m^ (46) 

r 

^DFT^^ = ^exp[z^-rl5^^^ (47) 

k 

Beware of the mapping between r and k, some routines use a weird wrapping! 

There are different ways of placing galaxies (or particle in your simulation) on a grid: 
Nearest grid point, Cloud in cell, trangular shaped cloud etc... For each of these remember{\) 
then to deconvolve the resulting P{k) for their effect. Note that 

4^(^)'^%--:.^^?- (48) 

and thus 

The discretization introduces several effects: 
The Nyquist frequency: k^y = is that of a mode which is sampled by 2 grid points. 
Higher frequencies cannot be properly sampled and give aliasing (spurious transfer of power) 
effects. You should always work at k < kNy There is also a minimum k (largest possible 
scale) that you finite box can test :kmin > 27r/L. This is one of the -many- reason why one 
needs ever larger N-body simulations... 

In addition DFT assume periodic boundary conditions, if you do not have periodic bound- 
ary conditions then this also introduces aliasing. 
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2. Window, selection function, masks etc 

Selection function: Galaxy surveys are usually magnitude limited, which means that 
as you look further away you start missing some galaxies. The selection function tells you 
the probability for a galaxy at a given distance (or redshift z) to enter the survey. It is a 
multiplicative effect along the line of sight in real space. 

Window or mask You can never observe a perfect (or even better infinite) squared box 
of the Universe and in CMB studies you can never have a perfect full sky map (we live in 
a galaxy...). The mask (sky cut in CMB jargon) is a function that usually takes values of 
or land is defined on the plane of the sky (i.e. it is constant along the same line of sight). 
The mask is also a real space multiplication effect. In addition sometimes in CMB studies 
different pixels may need to be weighted differently, and the mask is an extreme example of 
this where the weights are either or 1. Also this operation is a real space multiplication 
effect. 

Let's recall that a multiphcaton in real space (where W{x) denotes the effects of window 
and selection functions) 

(5*™^(f) — > 5°^'{x) = 5^'''^{x)W{x) (50) 
is a convolution in Fourier space: 

^true^^^ , d°^'{k) = 5*™"(fc) * W{k) (51) 

the sharper W{f) is the messier and delocahzed W{k) is. As a result it will couple different k- 
modes even if the underlying ones were not correlated! Discreteness While the dark matter 
distribution is almost a continuous one the galaxy distribution is discrete. We usually assume 
that the galaxy distribution is a sampling of the dark matter distribution. The discreteness 
effect give the galaxy distribution a Poisson contribution (also called shot noise contribution). 
Note that the Poisson contribution is non Gaussian: it is only in the limit of large number 
of objects (or of modes) that it approximates a Gaussian. Here it will suffice to say that 
as long as a galaxy number density is high enough (which will need to be quantified and 
checked for any practical application) and we have enough modes, we say that we will have 
a superposition of our random field (say the dark matter one characterized by its P{k)) plus 
a white noise contribution coming from the discreteness which amplitude depends on the 
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average number density of galaxies (and should go to zero as this go to infinity), and we 
treat this additional contribution as if it has the same statistical properties as the underlying 
density field (which is an approximation). What is the shot noise effect on the correlation 
properties? 



Following (1231 ) we recognize that our random field is now given by 

f{x) = n{x) =n[l + 6{x)] = ^ 5^(f - (52) 

i 

where n denotes average number of galaxies: n =< ~ ^i) >■ Then, as done when 

introducing the Poisson distribution, we divide the volume in infinitesimal volume elements 
6V so that their occupation can only be or 1. For each of these volumes the probability of 
getting a galaxy is 6P = p{x)6V, the probability of getting no galaxy is SP = 1 — p{x)6V 
and < rii >=< nf >= n6V. We then obtain a double stochastic process with one level of 
randomness coming from the underlying random field and one level coming from the Poisson 
sampling. The correlation function is obtained as: 

(E ^""(^"1 - rDS'^ir-^ - r-)) = n\l + ^12) + nd^'in - r^s) (53) 

thus 

< nin2 >= < S1S2 >'^] where < S1S2 >'^= ^{xu) + ^5^(rr-r"^) (54) 

and in Fourier space 

< Sk,Sk, {27rf (^P{k) + S'^ih + h) (55) 

This is not a complete surprise: the power spectrum of a superposition of two independent 
processes is the sum of the two power spectra.... 



3. pros and cons of £,{r) and P{k) 

Let us briefiy recap the pros and cons of working with power spectra or correlation 
functions. 

Power spectra 

Pros: Direct connection to theory. Modes are uncorrelated (in the ideal case). The average 
density ends up in P{k = 0) which is usually discarted, so no accurate knowledge of the 
mean density is needed. There is a clear distinction between linear and non-linear scales. 
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Smoothing is not a problem (just a multiplication) 

Cons: window and selection functions act as complicated convolutions, introducing mode 
coupling! (this is a serious issue) 
Correlation function 

Pros: no problem with window and selection function 

Cons: scales are correlated, covariance calculation a real challenge even in the ideal case. 
Need to know mean densities very well. No clear distinction between linear and non-linear 
scales. No direct correspondence to theory. 

4. ... and for CMB? 

If we can observe the full sky the the CMB temperature fluctuation field can be nicely 
expanded in spherical harmonics: 

e 

AT(n) = E E (^imYemin). (56) 
e>o m=-e 

where 

aem = / dnnAT{h)Y;^{h). (57) 

and thus 

< |Ofm|^ >= {(^emO-e'm') — ^U'Smm'Ci (58) 

Q is the angular power spectrum and 

Now what happens in the presence of real world effects such as a sky cut? Analogously 
to the real space case: 

aem = / dn,,AT{n)Win)Y;jn) (60) 

where W{n) is a position dependent weight that in particular is set to on the sky cut. 
As any CMB observation gets pixelized this is 

aem = ^pY: AT{p)W{p)Y;Jp) (61) 
p 

where p runs over the pixels and flp denotes the solid angle subtended by the pixel. 
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Clearly this can be a problem (this can be a nasty convolution), but let us initially ignore 
the problem and carry on. 

The pseudo-Q' s (Hivon et al 2002) (Q are defined as: 

1 ^ 



(2£ + 1) 



m=-i 



Clearly Ce 7^ Q but 



{Ce) = J2Gie'{Ce') (63) 



where <> denotes the ensamble average. 

We notice already two things: as expected the effect of the mask is to couple otherwhise 
uncorrelated modes. In large scale structure studies usually people stop here: convolve 
the theory with the various real world effects including the mask and compare that to the 
observed quantities. In CMB usually we go beyond this step and try to deconvolve the real 
world effects. 

First of all note that 

Ge.e. = E(24 + l)We, (|,^[,^'of (64) 

^3 



where 

We = IW'^-I' and W^m = f dQ^W {n)Y;^{h) (65) 

So if you are good enough to be able to invert G and you can say that < > is the 
you want then 

Ce = Y.Gie}Ce' (66) 
Not for all experiments it is viable (possible) to do this last step. 

In adddition to this, the instrument has other effects such as noise and a finite beam. 



5. Noise and beams 

Instrumental noise and the finite resolution of any experiment affect the measured C^. 
The effect of the noise is easily found: the instrumental noise is an independent random 
process with a Gaussian distribution superposed to the temperature field. In aim space 

„ , signal . noise 
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While (a^;^^^) — 0, in the power spectrum this gives rise to the so-called noise bias: 

^measured ^signal _|_ ^noise (6T) 

where Cf"'^'^ = i{2£ + 1) J^m k^m^'^P- As the expectation value of C*^"**^ is non zero, this is 
a biased estimator. 

Note that the noise bias disappears if one computes the so-called cross Ci obtained 
as a cross-correlation between different, uncorrelated, detectors (say detector a and b) as 
(^"m*^'"^"m*^''') — 0- is however not getting something for nothing: when one computes 
the covariance (or the associated error) for auto and for cross correlation Q (exercise!) the 
covariance is the same and includes the extra contribution of the noise. It is only that the 
cross-Q are unbiased estimators. 

Every experiment sees the CMB with a finite resolution given by the experimental beam 
(similar concept to the Point Spread Function for optical astronomy). The observed tem- 
perature field is smoothed on the beam scales. Smoothing is a convolution in real space: 

Ti^ J dQ!j'{h)b{\n - n'\) (68) 

where we have considered a symmetric beam for simplicity. The beam is often well ap- 
proximated by a Gaussian of a given Full Width at Half Maximum. Remember that 
(76 = QA2bFWHM. 

Thus in harmonic space the beam effect is a multiplication: 

(jmeasured ^ cf^^'^^^l (69) 

and in the presence of instrumental noise 

^measured ^^ky ^—(!^a'^ _|_ ^noise (70) 

Of course, one can always deconvolve for the effects of the beam to obtain an estimate of 

^measured ^ ^^^gg ^ possible tO Cf^: 

^measured' ^sky _|_ ^noise^Pa^ (71) 

That is why it is often said that the effective noise "blows up" at high I (small scales) 
and why it is important to know the beam(s) well. 
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Exercise: 
Exercise: 



What happens if you use cross-Q's? 

What happens to (7^easured' ^^^e beam is poorly reconstructed? 



Note that the signal to noise of a CMB map depends on the pixel size (by smoothing the 
map and making larger pixels the noise per pixel will decrease as ^^pix, ^pix being the new 
pixel solid angle), on the integration time api^ = s/y/t where s is the detector sensitivity 
and t the time spent on a given pixel and on the number of detectors api^ = s/\fM where 
M is the number f detectors. 

To compare maps of different beam sizes it is useful to have a noise measure that in 
independent of Vtpi^: w = (ali^Qpi^)'^ . 



Exercise: Compute the expression for (7^°*^*^ given: 
t = observing time 
s = detector sensitivity (in fxK/ ^/s) 
n = number of detectors 

= number of pixels 
fsky = fraction of the sky observed 



Assume uniform noise and observing time uniformly distributed. You may find (jlSl ) very 
useful. 



6. Aside: Higher orders correlations 

From what we have learned so far we can conclude that the power spectrum (or the 
correlation function) completely characterizes the statistical properties of the density field 
if it is Gaussian. But what if it is not? 

Higher order correlations are defined as: {6i...6m) where the deltas can be in real space 
giving the correlation function or in Fourier space giving power spectra. 

At this stage, it is useful to present here the Wick's theorem (or cumulant expansion 
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theorem). The correlation of order m can in general be written as sum of products of 
unreducible (connected) correlations of order £ ior i — l...m. For example for order 3 we 
obtain: 



(616263) f = 

{61) {62) (63) + 

{61) {6263) + (Scycterms) 

{616263) 



(72) 
(73) 
(74) 
(75) 



and for order 6 (but for a distribution of zero mean): 



< 6i..6e > 



< 61. ..6e >/= 

< 6162 >< 6364 >< 656Q > +..{15terms) 



< S1S2 >< 636^1656^ > +..{15terms) 



< 616263 >< 64656G > +...{10terms) 



(76) 
(77) 
(78) 
(79) 
(80) 



For computing covariances of power spectra, it is useful to be familiar with the above 
expansion of order 4. 

V. MORE ON LIKELIHOODS 

While the CMB temperature distribution is Gaussian (or very close to Gaussian) the 
distribution is not. At high i the Central Limit Theorem will ensure that the likelihood is 
well approximated by a Gaussian but at low i this is not the case. 



where T denotes a vector of the temperature map, Cj^ denotes the Ci given by a theoretical 
model (e.g. a cosmological parameters set), and Sij is the signal covariance: 




exp[-(T5-iT)/2] 



(81) 




(2£+l) 
An 



(82) 



and Pi denote the Legendre polynomials. 
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If we then expand T in spherical harmonics we obtain: 

£(T|Cf ) « 



exp[-l/2|a,^|VCf ] 



(83) 



Isotropy means that we can sum over m's thus: 



-21n£ = ^(2^+l) 



In 



C, 



th 



cf 



ata 



+ 



'c: 



data ^ 



Cf 



- 1 



(84) 



where Cf = Em |a€m|V(2^ + 1)- 



Exercise: show that for an experiment with (gaussian) noise the expression is the same but 
with the substitution C^'* — > Cj^ + A/^ with A/" denoting the power spectrum of the noise. 

Exercise: show that for a partial sky experiment (that covers a fraction of sky fsky you 
can approximately write: 

ln>C — > fsky'^nC (85) 
Hint: Think about how the number of independent modes scales with the sky area. 



As an aside... you could ask: "But what do I do with polarization data?". Well... if 
the aj^ are Gaussianly distributed also the af^ and af^ will be. So we can generalize the 
approach above using a vector (aJm'^fm'^fm) ■ Let us consider a full sky, ideal experiment. 
Start by writing down the covariance, follow the same steps as above and show that: 

( / r<BB\ /r<TTr<EE (r<TE\2\ 

-21n£ = V(2£ + 1) In fe- + In ' ' ~ ^ / J 



e 



r'TTr'EE I nTTnEE onTE/^TE r'BB 
r<TTr<EE (nTE\2 r'BB 



where Ce denotes Cf^ and Ce denotes Q'"*". 

It is easy to show that for a noisy experiment then Cf^ — > Cf^ + A/"/"^ where A/"^ 

denotes the noise power spectrum and X,Y = {T. E, B}. 

Exercise: generalize the above to partial sky coverage: for added complication take ^ 
ffky 7^ fsky (^^i^ often the case as the sky cut for polarization may be different from that 
of temperature( the foregrounds are different) and in general the cut (or the weighting) for 
B may need to be larger than that for E. 
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Following ( 1271 ) let us now expand in Taylor series Equation around its maximum by 



writing Ci = C\ {1 + e). For a single multipole 



'^2 ^3 



- 2 ln£, = {2^ + l)[e - ln(l + e)] ~ (2£ + 1) - - + 0{e^)^ . (87) 

We note that the Gaussian likelihood approximation is equivalent to the above expression 
truncated at e^: -2 In^causs,^ oc (2£ + 1)/2[(Q - Cf)/Cff ~ (2£ + l)eV2. 

Also widely used for CMB studies is the lognormal likelihood for the equal variance 
approximation (Bond et al 1998): approximation is 

In 



21n£; -(2^+^) 



'^2 ^3^ 



Thus our approximation of likelihood function is given by the form, 



1 2 

In £ = - In ^Gauss + 3 '^LN ' (89) 



where 



In^Gauss OC Y^^Cf - C,)Qu'{Cf - C,) , (90) 

and 

ln£LN = -l/2E(-2f - - ^^'), (91) 

where z^*^ = ln(C^'^ + A/"^), ^£ = ln(C£ + A/"^) and Q^/ is the local transformation of the 
curvature matrix Q to the lognormal variables z^, 

Q.,e = (Cf + Afi)Qw{Cj^ + Me). (92) 

The curvature matrix is the inverse of the covariance matrix evaluated at the maximum 
likelihood. However we do not want to adopt the "equal variance approximation", so Q for 
us will be in inverse of the covariance matrix. 

The elements of the covariance matrix, even for a ideal full sky experiment can be written 

as: 



^« = ^WTi ^^^^ 
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In the absence of noise the covariance is non zero: this is the cosmic variance. 



Note that for the latest WMAP release (j22 



26l)at low i the likelihood is computed directly 



from the maps rh. The standard likelihood is given by 



L{Th\S)drh 



exp 



dm 



(94) 



\S + N\y^ (27r)3"p/2' 

where rh is the data vector containing the temperature map, T, as well as the polarization 
maps, Q, and U, Up is the number of pixels of each map, and 5* and are the signal and 
noise covariance matrix {Sup x Sup), respectively. As the temperature data are completely 
dominated by the signal at such low multipoles, noise in temperature may be ignored. This 
simplifies the form of likelihood as 



L{rh\S)dfh 



exp 


-^m{Sp + Np) 




dm 


exp(-if*5^if) df 




\Sp + Np\y^ 




(27r)"p 


IStI^/^ (27r)"p/2 



(95) 



where St is the temperature signal matrix {up x Up), the new polarization data vector, 
m = {Qp, Up) and Sp is the signal matrix for the new polarization vector with the size of 

2 Tip X 2 Tip. 

At the time of writing, in CMB parameter estimates for i < 2000, the likelihood calcu- 
lation is the bottleneck of the analysis. 



VI. MONTE CARLO METHODS 
A. Monte Carlo error estimation 



Let's go back to the issue of parameter estimation and error calculation. Here is the 
conceptual interpretation of what it means that en experiment measures some parameters 
(say cosmological parameters). There is some underlying true set of parameters dtme that 
are only known to Mother Nature but not to the experimenter. There true parameters 
are statistically realized in the observable universe and random measurement errors are 
then included when the observable universe gets measured. This "realization" gives the 
measured data Vq. Only Vq is accessible to the observer (you). Then you go and do what 
you have to do to estimate the parameters and their errors (chi-square, likelihood, etc..) 
and get do. Note that Vq is not a unique realization of the true model given by dtme '■ 
there could be infinitely many other realizations as hypothetical data sets, which could have 
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been the measured one: Vi,V2, ... each of them with a shghtly different fitted parameters 
di,d2---- (So is one parameter set drawn from this distribution. The hypotetical ensamble 
of universes described by di is called ensamble, and one expects that the expectation value 
(di) = dtrue- If we knew the distribution of di — dtme we would know everything we need 
about the uncertainties in our measurement d^. The goal is to infer the distribution of 

di - dtrue without kuowiug dtrue- 

Here's what we do: we say that hopefully ao is not too wrong and we consider a fictitious 
world where was the true one. So it would not be such a big mistake to take the 
probability distribution of ao ~ di to be that of dtrue — di. In many cases we know how to 
simulate d^ — di and so we can simulate many synthetic realization of "worlds" where is 
the true underlying model. Then mimic the observation process of these fictitious Universes 
replicating all the observational errors and effects and from each of these fictitious universe 
estimate the parameters. Simulate enough of them and from af — cJo you will be able to 
map the desired multi-dimensional probability distribution. 

With the advent of fast computers this technique has become increasingly widespread. 
As long as you believe you know the underlying distribution and that you believe you can 
mimic the observation replicating all the observational effects this technique is extremely 
powerful and, I would say, indispensable. 



B. Monte Carlo Markov Chains 



When dealing with high dimensional likelihoods (i.e. many parameters) the process of 
mapping the likelihood (or the posterior) surface can become very expensive. For example for 
CMB studies the modles considered have from 6 to 11+ parameters. Every model evaluation 
even with a fact code such as CAMB can take up to minutes per iteration. A grid-based 
likelihood analysis would require prohibitive amounts of CPU time. For example, a coarse 
grid (~ 20 grid points per dimension) with six parameters requires ~ 6.4 x 10^ evaluations 
of the power spectra. At 1.6 seconds per evaluation, the calculation would take ~ 1200 
days. Christensen & Meyer (2000) (jj) proposed using Markov Chain Monte Carlo (MCMC) 
to investigate the likelihood space. This approach has become the standard tool for CMB 
analyses. MCMC is a method to simulate posterior distributions. In particular one simulates 
sampling the posterior distribution 'P(alx), of a set of parameters a given event x, obtained 
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via Bayes' Theorem 

^/ , X V(x\a)V(a) 
^ ' ^ JP{x\a)V{a)da ^ ' 

where V{x\a) is the hkehhood of event x given the model parameters a and 'P(a) is the prior 

probabihty density; a denotes a set of cosmological parameters (e.g., for the standard, flat 

ACDM model these could be, the cold-dark matter density parameter f2c, the baryon density 

parameter the spectral slope n^, the Hubble constant -in units of 100 km s~^ Mpc~^)- 

/i, the optical depth r and the power spectrum amplitude A), and event x will be the set 

of observed Ci. The MCMC generates random draws (i.e. simulations) from the posterior 

distribution that are a "fair" sample of the likelihood surface. From this sample, we can 

estimate all of the quantities of interest about the posterior distribution (mean, variance, 

confidence levels). The MCMC method scales approximately linearly with the number of 

parameters, thus allowing us to perform likelihood analysis in a reasonable amount of time. 

A properly derived and implemented MCMC draws from the joint posterior density 
V{a\x) once it has converged to the stationary distribution. The primary consideration in 
implementing MCMC is determining when the chain has converged. After an initial "hum- 
in" period, all further samples can be thought of as coming from the stationary distribution. 
In other words the chain has no dependence on the starting location. 

Another fundamental problem of inference from Markov chains is that there are always 
areas of the target distribution that have not been covered by a finite chain. If the MCMC 
is run for a very long time, the ergodicity of the Markov chain guarantees that eventually 
the chain will cover all the target distribution, but in the short term the simulations cannot 
tell us about areas where they have not been. It is thus crucial that the chain achieves 
good "mixing". If the Markov chain does not move rapidly throughout the support of the 
target distribution because of poor mixing, it might take a prohibitive amount of time for 
the chain to fully explore the likelihood surface. Thus it is important to have a convergence 
criterion and a mixing diagnostic. Plots of the sampled MCMC parameters or likelihood 
values versus iteration number are commonly used to provide such criteria (left panel of 
Figure [3]). However, samples from a chain are typically serially correlated; very high auto- 
correlation leads to little movement of the chain and thus makes the chain to "appear" to 

n 

have converged. For a more detailed discussion see (110). Using a MCMC that has not fully 
explored the likelihood surface for determining cosmological parameters will yield wrong 
results. See right panel of Figure [3]). 
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C. Markov Chains in Practice 

Here are the necessary steps to run a simple MCMC for the CMB temperature power 
spectrum. It is straightforward to generahze these instructions to include the temperature- 
polarization power spectrum and other datasets. The MCMC is essentially a random walk in 
parameter space, where the probability of being at any position in the space is proportional 
to the posterior probability. 

1) Start with a set of cosmological parameters {ai}, compute the Cj and the hkelihood 

Ci = C{Cj^^^\C£). 

2) Take a random step in parameter space to obtain a new set of cosmological parameters 
{a2}. The probability distribution of the step is taken to be Gaussian in each direction 
i with r.m.s given by (Xj. We will refer below to o"j as the "step size". The choice of 
the step size is important to optimize the chain efficiency (see discussion below) 

3) Compute the C^*^ for the new set of cosmological parameters and their hkelihood £2- 

4. a) If £2/^1 > 1, "take the step" i.e. save the new set of cosmological parameters {a2} 
as part of the chain, then go to step 2 after the substitution {ai} — > {0^2} ■ 

4.b) If >C2/>Ci < 1, draw a random number x from a uniform distribution from to 1. If 
X > jC,2/jCi "do not take the step", i.e. save the parameter set {ai} as part of the 
chain and return to step 2. li x < JC2/C1, " take the step", i.e. do as in 4. a). 

5) For each cosmological model run four chains starting at randomly chosen, well- 
separated points in parameter space. When the convergence criterion is satisfied and 
the chains have enough points to provide reasonable samples from the a posteriori 
distributions (i.e. enough points to be able to reconstruct the 1- and 2-a levels of the 
marginalized likelihood for all the parameters) stop the chains. 

It is clear that the MCMC approach is easily generalized to compute the joint likefihood of 
WMAP data with other datasets. 
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D. Improving MCMC Efficiency 



MCMC efficiency can be seriously compromised if there are degeneracies among parame- 
ters. The typical example is the degeneracy between Q^n and H for a flat cosmology or that 
between fi^ and for a non-flat case (see e.g. reference (H)). 

The Markov chain efficiency can be improved in different ways. Here we report the 
simplest way. 

Repar amet er izat ion 

We describe below the method we use to ensure convergence and good mixing. Degen- 
eracies and poor parameter choices slow the rate of convergence and mixing of the Markov 
Chain. There is one near-exact degeneracy (the geometric degeneracy) and several ap- 
proximate degeneracies in the parameters describing the CMB power spectrum Bond et al 
(1994) (jl|), Efstathiou& Bond(1984) (^. The numerical effects of these degeneracies are 
reduced by finding a combination of cosmological parameters (e.g., Qc, ^b, h, etc.) that 
have essentially orthogonal effects on the angular power spectrum. The use of such pa- 
rameter combinations removes or reduces degeneracies in the MCMC and hence speeds up 
convergence and improves mixing, because the chain does not have t o sp end time exploring 
degeneracy directions. Kosowsky, Milosavljevic & Jimenez (2002) p9!) and Jimenez et al 
(2003) (jl6l ) introduced a set of reparameterizations to do just this. In addition, these new 
parameters reflect the underlying physical effects determining the form of the CMB power 
spectrum (we will refer to these as physical parameters). This leads to particularly intuitive 
and transparent parameter dependencies of the CMB power spectrum. 

For the 6 parameters LCDM model these "normal" or "physical" parameters are: the 
physical energy densities of cold dark matter, ujc = ^ch"^, and baryons, Ub = flbh"^, the 
characteristic angular scale of the acoustic peaks, 

Oa = (97) 

where a^ec is the scale factor at decoupling, 

rsiadec) = X (98) 

"■dec dx 



1 + ^) ((1 - + + n^X + nrad) ^'^ 
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is the sound horizon at decouphng, and 



dAM = ^^fk^ (99) 
^0 J\Q - 11 



where 

dx 



r(adec) = \^~M (100) 

[(1 - fi)x2 + + fi^a; + a,,]-'/' 

and S^ir) as usual coincides with the argument if the curvature k is 0, is a sin function 
for n > 1 and a sinh function otherwhise. Here Hq denotes the Hubble constant and c is 
the speed of light, Qm = + ^b, denotes the dark energy density parameters, w is the 
equation of state of the dark energy component, Q = Qm + and the radiation density 
parameter Q^ad = ^-y + ^-y, are the the photon and neutrino density parameters 
respectively. For reionization sometmes the parameter Z = exp(— 2r) is used, where r 
denotes the optical depth to the last scattering surface (not the decoupling surface). 

These reparameterizations are useful because the degeneracies are non-linear, that is they 
are not well described by ellipses in parameter space. For degeneracies that are well approx- 
imated by ellipses in parameter space it is possible to find the best reparameterization auto- 
matically. This is what the code CosmoMC (5;!^ (see tutorials) does. To be more precise 
it computes the parameters covariance matrix from which the axes of the mult i- dimensional 
degeneracy ellipse can be found. Then it performs a rotation and re-scaling of the coordi- 
nates (i.e. the parameters) to transform the degeneracy ellipse in an azimutally symmetric 
contour. See discussion at http://cosmologist.info/notes/CosmoMC.pdf for more in- 
formation. This technique can improve the MCMC efficiency up to a factor of order 10. 

Step size optimization The choice of the step size in the Markov Chain is crucial 
to improve the chain efficiency and speed up convergence. If the step size is too big, the 
acceptance rate will be very small; if the step size is too small the acceptance rate will be 
high but the chain will exhibit poor mixing. Both situations will lead to slow convergence. 



E. Convergence and Mixing 

Before we even start this section: thou shall always use a convergence and mixing 
criterion when running MCMC's. 

Let's illustrate here the method proposed by Gelman & Rubin 1992 as an example. 
They advocate comparing several sequences drawn from different starting points and check- 
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ing to see that they are indistinguishable. This method not only tests convergence but can 
also diagnose poor mixing. Let us consider M chains starting at well-separated points in 
parameter space; each has 2A'" elements, of which we consider only the last N: {yj} where 
i = 1,..,N and j = 1,..,M, i.e. y denotes a chain element (a point in parameter space) 
the index i runs over the elements in a chain the index j runs over the different chains. We 
define the mean of the chain 

y'-^Y^yi, (101) 



and the mean of the distribution 



N . , 
1=1 



NM 



tj = l 



We then define the variance between chains as 

M 



Bn = Y.{y' - yf^ (103) 



and the variance within a chain as 

w ' 



M(N - 1) ^ 

The quantity 



w 



;v (105) 



is the ratio of two estimates of the variance in the target distribution: the numerator is an 
estimate of the variance that is unbiased if the distribution is stationary, but is otherwise an 
overestimate. The denominator is an underestimate of the variance of the target distribution 
if the individual sequences did not have time to converge. 

The convergence of the Markov chain is then monitored by recording the quantity R for 
all the parameters and running the simulations until the values for R are always < 1.03. 
Needless to say that the cosmomc package offers several convergence and mixing diagnostic 
tools as part of the "getdist" routine. 



Question: how does the MCMC sample the prior if all one actually computes is the 
likelihood? 
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F. MCMC Output Analysis 

Now that you have your multiple chains and the convergence criterium says they are 
converged whot do you do? First discard burn in and merge the chains. Since the MCMC 
passes objective tests for convergence and mixing, the density of points in parameter space 
is proportional to the posterior probability of the parameters. (Note that cosmomc saves 
repeated steps as the same entry in the file but with a weight equal to the repetitions: the 
MCMC gives to each point in parameter space a "weight" proportional to the number of steps 
the chain has spent at that particular location.). The marginalized distribution is obtained 
by projecting the MCMC points. This is a great advantage compared to the grid-based 
approach where multi-dimensional integrals would have to be performed. The MCMC basi- 
cally performs a Monte Carlo integration, the density of points in the n-dimensional space 
is proportional to the posterior, and best fit parameters and multi-dimensional confidence 
levels can be found as illustrated in the last class. 

Note that the global maximum likelihood value for the parameters does not necessarily 
coincide with the expectation value of their marginalized distribution if the likelihood surface 
is not a multi-variate Gaussian. 

A virtue of the MCMC method is that the addition of extra data sets in the joint analysis 
can efficiently be done with minimal computational effort from the MCMC output if the 
inclusion of extra data set does not require the introduction of extra parameters or does not 
drive the parameters significantly away from the current best fit. If the likelihood surface for 
a subset of parameters from an external (independent) data set is known, or if a prior needs 
to be added a posteriori, the joint posterior surface can be obtained by multiplying the new 
probability distribution with the posterior distribution of the MCMC output. To be more 
precise: as the density of point (i.e. the weight) is directly proportional to the posterior, 
then this is achieved by multiplying the weight by the new probability distribution. 

The cosmomc package already includes this facility. 

VII. FISHER MATRIX 

What if you wanted to forecast how well a future experiment can do? There is the 
expensive but more accurate way and the cheap and quick way (but often less accurate) . The 
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expensive way is in the same spirit of Monte-Carlos simulations discussed earlier: simulate 
the observations and estimate the parameters as you would do on the data. However often 
you want to have a much quicker way to forecasts parameters uncertainties, especially if you 
need to quickly compare many different experimental designs. This technique is the Fisher 
matrix approach. 



VIII. FISHER MATRIX 



The question of how accurately one can measure model parameters from a given data set 
(without simulating the data set) was answered more than 70 years ago by Fisher (1935) 
(j9|). Suppose your data set is given by m real numbers Xi...Xm arranged in a vector (they 
can be CMB map pixels, P{k) of galaxies etc..) x is a random variable with probability 
distribution which depends in some way on the vector of model parameters a. The Fisher 
information matrix id defined as: 

where L = —InC. In practice you will choose a fiducial model and compute the above at 
the fiducial model. In the one parameter case let's note that if the Likelihood is Gaussian 
then L = l/2{a — ao)/a'^ where ao is the value that maximizes the likelihood and sigma is 
the error on the parameter a. Thus the second derivative wrt a of L is 1/cr^ as long as aa 
does not depend on a. In general we can expand L in Taylor series around its maximum 
(or the fiducial model). There by definition the first derivative of L wrt the parmeters is 0. 

When 2AL = 1, which in the case when we can identify it with Ax^ with corresponds to 



68% (or 1-0") then 1/ ^(PL/da^ is the 1 sigma displacement of a from a^. 

The generalization to many variables is beyond our scope here (see Kendall & Stuard 



1977 (1l7l )) let's just say that an estimate of the covariance for the parameters is given by: 

> (F-^)., (108) 
From here it follows that if all other parameters are kept fixed 

(109) 

V ^ a 
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(i.e. the reciprocal of the square root of the diagonal element ii of the Fisher matrix.) 

But if all other parameters are estimated from the data as well then the marginalized 
error is 

= (F-^)f (110) 

(i.e. square root of the element ii of the inverse of the Fisher matrix -perform a matrix 
inversion here!) 

In general, say that you have 5 parameters and that you want to plot the joint 2D 
contours for parameters 2 and 4 marginalized over all other parameters 1,3,5. Then you 
invert Fij, take the minor 22, 24, 42, 24 and invert it back. The resulting matrix, let's call it 
Q, describes a Gaussian 2D likelihood surface in the parameters 2 and 4 or, in other words, 
the chisquare surface for parameters 2,4 - marginalized over all other parameters- can be 
described by the equation = J2kq{(^k — Oik'^^'^^'^^)Qkq{oiq — a^*'^"'^*"')- 

From this equation, getting the errors corresponds to finding the quadraitic equation 
solution = Ax^. For correspondence between Ax^ and confidence region see the earlier 
discussion. If you want ti make plots, the equation for the elliptical boundary for the joint 
confidence region in the sub-space of parameters of interest is: A = 5dQ~^5a. 

Note that in many applications the likelihood for the data is assumed to be Gaussian 
and the data-covariance is assumed not to depend on the parameters. For example for the 
application to CMB the fisher matrix is often computed as: 

_ (2£ + 1) nm 

This approximation is good in the high i, noise dominated regime. In this regime the 
central limit theorem ensures a Gaussian likelihood and the cosmic variance contribution 
to the covariance is negligible and thus the covariance does not depend on the cosmological 
parameters. However at low i in the cosmic variance dominated regime, this approximation 
over-estimates the errors (this is one of the few cases where the Fisher matrix approach 
over-estimates errors) by about a factor of 2. In this case it is therefore preferable to go 
for the more numerically intensive option of computing eq. 11061 with the exact form for the 
likelihood Eq[HSl 

Before we conclude we report here a useful identity: 

2L = lndet{C) + (x - y)^Cr^J\x - y)J = Tr[lnC + C-^D] (112) 
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where C stands for the covariance matrix of the data, repeated indices are summed over, 
X denotes the data and y the model fitting the data and D is the data matrix defined as 
{x — y){x — y)^. We have used the identity In det{C) ~ Tr In C. 

IX. CONCLUSIONS 

Wether you are a theorist or an experimentahst in cosmology, these days you cannot ignore 
the fact that to make the most of your data, statistical techniques need to be employed, and 
used correctly. An incorrect treatment of the data will lead to nonsensical results. A given 
data set can revel a lot about the universe, but there will always things that are beyond the 
statistical power of the data set. It is crucial to recognize it. I hope I have given you the 
basis to be able to learn this by yourself. When in doubt always remember: treat your data 
with respect. 
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X. QUESTIONS 

Here are some of the questions (and their answers...) 

Q: What about the forgotten V{D) in the Bayes theorem? 
A: V{D) becomes important when one wants to compare different cosmological models (not 
just different parameter values within the same cosmological model). There is a somewhat 
extensive literature in cosmology alone on this: "Baycsian evidence" is used to do this "model 
selection" . Bayesian evidence calculation can be, in many cases, numerically intensive. 

Q: What do you do when the likelihood surface is very complicated for example is multi- 
peaked? 

A: Luckily enough in CMB analysis this is almost never the case (exceptions being, for exam- 
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pie, the cases where sharp oscillations are present in the Cj^ as it happens in transplanckian 
models.) In other contexts this case is more frequent. In these cases, when additional peaks 
can't be suppressed by a motivated prior, there are several options: if m is the number 
of parameters, report the confidence levels by considering only regions with m-dimensional 
likelihoods above a threshold before marginalizing (thus local maxima will be included in 
the confidence contours if significant enough) or b) simply marginalize as described here, but 
expect that the marginalized peaks will not coincide with the m-dimensional peaks. The 
most challenging issue when the likelihood surface is complicated is having the MCMC to 
fully explore the surface. Techniques such as Hamiltonian Monte Carlo are very powerful in 
these cases see e.g. ( Illl : Il2l ) 

Q: The Fisher matrix approach is very interesting, is there anything I should look out 
for? 

A: The Fisher matrix approach assumes that the parameters log-likelihood surface can be 
quadratically approximated around the maximum. This may or may not be a good ap- 
proximation. It is always a good idea to check at least in a reference case wether this 
approach significantly underestimated the errors. In many Fisher matrix implementation 
the likelihood for the data is assumed to be Gaussian and the data-covariance is assumed 
not to depend on the parameters. While this approximation greatly simplifies the calcula- 
tions (Exercise: show why this is the case), it may significantly mis-estimate the size of 
the errors. In addition, the Fisher matrix calculation often requires numerical evaluation of 
second derivatives. Numerical derivatives always need a lot of care and attention: you have 
been warned. 

Q: Does cosmomc use the sampling described here? 
A: The recipe reported here is the so called Metropolis-Hasting algorithm. Cosmomc offers 
also other sampling algorithms: slice sampling, or a split between "slow" and "fast" parame- 
ters or the learn-propose option where chains automatically adjust the proposal distribution 
by repeatedly computing the parameters covariance matrix on the fly. These techniques can 
greatly improve the MCMC "efficiency". See for example 
http://cosmologist.info/notes/CosmoMC.pdf or (lid ) for more details. 
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XI. TUTORIALS 



Two tutorial sessions were organized as part of the Winter School. You may want to try 
to repeat the same steps. 

The goals of the first tutorial were: 



Download and install Healpix (|7|; Il4l ). make sure you have a fortran 90 compiler 
installed and possibly also IDL installed 



The LAMBDA site (1201 ) contains a lot of extremely useful information: browse the 
page. 

Find the on-line calculator for the theory Ce given some cosmological parameters, 
generate a set of Ci and save them as a "fits" file. 

Browse the help pages of Healpix and find out what it does. 

In particular the routine " symfast" enables one to generate a map from a set of theory 
Ce. 

The routine "anafast" enables one to compute Ce from a map. 

Using " symfast" , generate two maps with two different random number generators for 
the Ce you generated above. Select a beam of FWHM of, say, 30', but for now do 
not imose any galaxy cut and do not add noise. Compare the maps. To do that use 
"mollview". Why are they different? 



45 



Using "anafast" compute the Ci for both maps. 

Using your favorite plotting routine (the IDL part of healpix offers utihties to do that) 
plot the orginal Ce and the two realizations. Why do they differ? 

Deconvolve for the beam, and make sure you are plotting the correct units, the correct 
factors of i{i + 1) etc. Why do the power spectra still differ? 

Try to compute the Ci for imax = 3000 and for a non flat model with 
camb. It takes much longer than for a simple, fiat LCDM model. A 



code like CMBwarp ( ll6l ) offers a shortcut. try it by downloading it at 
http://www.astro.princeton.edu/~raulj/CMBwarp/index.htnil. Keep in mind 
however that this offers a fast fitting routine for a fine grid of pre-computed Ci, it is 
not a Boltzmann code and thus is valid only within the quoted limits! 

The goals of the second tutorial were: 

• Download and install the CosmoMC (5) package (and read the instructions). 

• Remember that WMAP data and likelihood need to be dowloaded separately from the 
LAMBDA site. Do this following the instructions. 

• Take a look at the params.ini file. In here you set up the MCMC 

• Set a chain to run. 

• get familar with the getdist program and distparams.ini file. This program checks 
convergence for you and compute basic parameter estimation from converged chains. 

• download from LAMBDA a set of chains, the suggested one was for the fiat, 
quintessence model for WMAP data only. 

• plot the marginalized probability distribution for the parameter w. 

• the challenge: Apply now a Hubble constant prior to this chain, take the HST key 
project (isl) constraint oi Hq = 72 ± 8 km/s/Mpc and assume it has a Gaussian 
distribution. This procedure is called "importance sampling". 
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FIG. 1 Example of linear fit to the data points Di in 2D. 
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Periodic boundary conditions 




Smooth transition 




FIG. 2 The importance of periodic boundary conditions 
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FIG. 3 Unconverged Markov chains. The left panel shows a trace plot of the likelihood values 
versus iteration number for one MCMC (these are the first 3000 steps from run). Note the burn-in 
for the first 100 steps. In the right panel, red dots are points of the chain in the (n, A) plane 
after discarding the burn-in. Green dots are from another MCMC for the same data-set and the 
same model. It is clear that, although the trace plot may appear to indicate that the chain has 
converged, it has not fully explored the likelihood surface. Using either of these two chains at this 



stage wil 



from ref (|27l ) 



_give incorrect results for the best fit cosmological parameters and their errors. Figure 
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TABLE I Ax^ for the conventionals 1, 2,and 3 — cr as a function of the number of parameters for 



the joint confidence levels. 
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1 
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1-a 


68.3% 


1.00 


2.30 


3.53 




90% 


2.71 


4.61 


6.25 


2-a 


95.4% 


4.00 


6.17 


8.02 


3-a 


99.73%: 


9.00 


11.8 


14.2 
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